Emergence of skew distributions in controlled growth processes 
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Starting from a master equation, we derive the evolution equation for the size distribution of 
elements in an evolving system, where each element can grow, divide into two, and produce new 
elements. We then probe general solutions of the evolution equation, to obtain such skew distri- 
butions as power-law, log-normal, and Weibull distributions, depending on the growth or division 
and production. Specifically, repeated production of elements of uniform size leads to power-law 
distributions, whereas production of elements with the size distributed according to the current 
distribution as well as no production of new elements results in log-normal distributions. Finally, 
division into two, or binary fission, bears Weibull distributions. Numerical simulations are also 
carried out, confirming the validity of the obtained solutions. 

PACS numbers: 05.40.-a, 89.75.Fb, 05.65.+b 



I. INTRODUCTION 

The central limit theorem is one of the most impor- 
tant theorems in probability theory. It states that the 
sum of a large number of independent random variables 
has a Gaussian distribution, which is symmetric. In re- 
ality, however, there exist many phenomena to which the 
central limit theorem is not applicable. In most of those 
cases, asymmetric skew distributions arise, characterized 
by long tails on one side; prototypic examples include 
power-law, log-normal, and Weibull distributions. In 
contrast to the Gaussian distribution, of which the origin 
is well known, there still lacks appropriate understanding 
of the origin of skew distributions. 

The power-law distribution is observed in various phys- 
ical, biological, economical, or social systems and there 
are extensive studies attempting to understand its prop- 
erties P-- 3]; inter alia, scale invariance is the most impor- 
tant property. While Zipf's law Q and Pareto's law [f| 
provide well-known examples of the power-law distribu- 
tion, the Yule process @, Gibrat's law [tJ, or preferential 
attachment @ has been proposed as the mechanism for 
the power-law distribution. Related to the power-law dis- 
tribution, the log-normal distribution has been used as 
an alternative of the former [914IT1]. Unlike the Gaussian 
distribution, the log-normal distribution emerges from 
multiplication (rather than addition) of a large number 
of independent random variables [12| ]. Thus geometric 
Brownian motion or the Black-Scholes equation is related 
to the log-normal distribution 13]. Interrelations be- 
tween the power-law distribution and the log-normal one 
have drawn a great deal of attention in literature and the 
concept of the double Pareto distribution has also been 
proposed 0, EH] • The Weibull distribution [l6[ is also 
one of the widely applie d distribution functions in a va- 
riety of systems (lU Il7i - ll9| . 
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distribution has been attributed to diverse origins on the 
case-by-case basis, including the extreme value statistics 
[IS HIl' fractal cracking fl 7| . Fickian diffusion in the frac- 
tal and non-fractal spacespa 23] , and the branching pro- 
cess [24| . Interrelations between the Weibull distribution 
and the log-normal one have also been studied and the 
similarity between them has been elucidated 25j. There 
are also works which analyzed a variety of time series 
data, e.g., from respiratory dynamics, electroencephalo- 
gram (EEG), heartbeats, and stock price records, and 
probed the asymmetry of distributions 26 29]. 

Even though some relations between these skew distri- 
butions are available, there still lacks theoretical under- 
standing of how the skew distributions appear in evolving 
systems. Recently, a general framework has been pro- 
posed, based on the master equation describing evolv- 
ing systems, to elucidate the emergence of skew distri- 
butions [3fj| . Following this approach, we begin with the 
same master equation, to obtain the time evolution equa- 
tion for the size distribution in evolving systems, which 
provides refinements and extensions. It is found that 
such skew distributions indeed emerge, depending on the 
growth, division, and production. In particular, the de- 
tailed conditions for power-law, log-normal, and Weibull 
distributions are clarified. In addition, numerical simu- 
lations are also carried out extensively; results are shown 
to agree well with the analytical ones. 

There are five sections in this paper: In Sec. II we 
introduce the master equation approach to the prolifera- 
tion process, including growth, division, and production. 
Section III is devoted to solving the evolution equation 
for various cases. Such skew distributions as power-law, 
log-normal, and Weibull distributions are obtained in the 
systems with uniform production, self-production, and 
binary fission. Section IV presents the results of numer- 
ical simulations, which confirm the analytical results in 
Sec. III. Finally, a summary together with a discussion 
of the results is given in Sec. V. 
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II. MASTER EQUATION 

We consider a system of N elements, the ith of which 
is characterized by the size x^ (i = 1, . . . , N). The con- 



figuration of the system is specified by the sizes of all 
elements, {x 1} . . . , x N }. The probability P(x 1 , . . . , x N ; t) 
for the system to be in configuration {x 1: . . . , x N } at time 
t is governed by the master equation 



d N r 

—P(x 1 ,...,x N ;t) = Y1J dx 'i ~* x t )P(x 1 ,...,x' t ,... ,x N ;t) - u{ Xl -+ x' l )P(x 1 , . . . , x N ;t)} , (1) 



where uj(xi -+ x^) is the transition rate for the ith ele- 
ment to change its size from xi to x\. Here we consider 
the case that the size changes by the amount propor- 
tional to the current size. The transition rate then takes 
the form 

uj(x -+ x') = X6(x' - (l+b)x) (2) 

with the (mean) growth rate A and the growth factor b 
30]. 

We are interested in the size distribution of the system 
f(x, t), related to the probability P(x 1 , . . . , x N ; t) via: 

1 f N 

f( x > t ) = ly / d N x ^2§(xi -x)P(x 1 ,...,x N ;t), (3) 

where J d N x = J °° dx 1 ■ ■ ■ J °° dx N . The evolution equa- 
tion for the size distribution f{x,t) can be obtained 



straightforwardly from Eq. (fTJ). In the case that the to- 
tal number of elements remains constant, the evolution 
equation for f(x,t) reads 

»/(,,,) —VfcO + j^jfj..). (4) 

Note the scale invariance of Eq. (Q} : Changing the scale 
according to x — >• cx and multiplying Eq. ((4]) by c for 
normalization, we obtain Eq. (|4]) in the form 

n c I cx \ 

— cf(cx, t) = -Xcf(cx, t) + A -/ — -, f ] . (5) 

dt Jy ' M ' ' 1 + 6 \l+b J y ' 

This manifests that for f(x, t) satisfying Eq. Q, cf(cx, t) 
also satisfies Eq. (j4]), thus providing a solution. 

On the other hand, when the total number of elements 
varies with time, i.e., N — N(t), it is not allowed to inter- 
change the order of d/dt and J d N x. To circumvent this 
difficulty, we deal with the time derivative of N(i)f(x, t): 



^-[N(t)f(x,t)] = lim i 
dV w v Jl At^o At 



. N+AN . N 

J d N+AN x S(x i -x)P(x 1 ,...,x N+AN ;t+At) -J d N xJ25(xi-x)P(x u ...,x N ;t) 



...,x N ; t+At) - P(x t , ...,x N ;t)} 

N+AN 



d x N+i ' ' ' dx N+AN 5{x i — x)P(x N+1 , . . . , x N+AN ; t+At) 

i=N+l 



(6) 



where it has been noted that J dx N+1 ■ ■ ■ dx N+AN P(x 1 , . . . ,x N+AN ;t+At) = P(x 1 , . . . , x N ; t+At) and 
J d N x P(x 1 , . . . , x N+AN ; t+At) — P(x N+1 , . . . , x N+AN ; t+At). 

I 

The first term on the right-hand side (the second line) and takes the form: 
of Eq. ^ contains the information only on the existing 
elements (up to N) and is exactly the same as the case 
for fixed N, thus resulting in 



(1st term) = N(t) 



-Xf(x,t) + 7T^/ ( 7T5>* 



(7) 



(2nd term) = N(t)g(x, t) = rN(t)g(x, t), (8) 



The second term (the third line) of Eq. ([6]) is involved 
with the information on the newly produced elements 



where g(x,t), define to be 



3 



j(x,t) 



1 

~AN 



dx 



N+l ' ' 



dx 



N+AN 



N+AN 
i=N + l 



x)P(x 



N+V ■ 



, X N+AN 



(9) 



r 



corresponds to the size distribution of newly produced 
elements at time t and r is the mean production rate. 
Namely, it is assumed that each element tends to pro- 
duce a new one with rate r, which in turn leads the total 
number of elements to increase in proportion to the cur- 
rent number: N = rN . Here care needs to be given to 
the limit At — > since AN is discrete, taking only integer 
values. The expression AN = rNAt is valid as long as 
NAt is finite, so that AN remains discrete. Accordingly, 
At can be reduced down to (rN) -1 , which corresponds 
to the time interval required to produce a new element 
(AN = 1). This vanishes as N grows large, thus allowing 
limAt^o AN/ At = N, etc. With the help of Eqs. ©-© 
and 

|/(M) = ^|[A7(M)]-^/0M), (10) 

the time evolution equation for the size distribution fi- 
nally obtains the form: 

(11) 



III. ASYMPTOTIC SOLUTIONS 

In this section we probe analytically the solutions of 
the evolution equations (j4]) and (|TTj) . In the case of 
no production (r = 0), the time evolution of f(x,t) is 
described by Eq. (U). Changing the variable from x 
to In a: = X, we write the distribution function for X: 
F(X,t) = e x f(e x ,t). Correspondingly, Eq. ((4]) becomes 



0_ 

di 



F(X,t) = -XF(X,t) + XF(X-a,t), (12) 



with a = ln(l+6). Given the initial distribution f(x, 0) = 
S(x - x ) or equivalently, F(X, 0) = 5(X - X ), F(X, t) 
vanishes unless X = X n = X +an for any integer n (note 
that this corresponds simply to x n — xoe an = (l+b) n Xo). 
With the definition F n (t) = F(X +an, t), Eq. JTJJ reads 



F n (t) = -XF n (t) + XF n _ 1 (t), 
which bears the Poisson distribution 



1 



F n (t) = - (Xt) , 



n -\t 



(13) 



(14) 



for n > and F n (t) = otherwise [31(. 

This result can be also recovered by considering the 
Fourier transform F(k,t) of Eq. (p~3|) . which satisfies 



d t F(k,t) = X{e 



-ika 



l)F(k,i). 



(15) 



The solution F(k, t) = exp[-A(l - er lka )t]F(k, 0) of Eq. 
([TB]) leads directly to 

F(X,t)= J dX'G F (X,X';t)F(X',0) 

= G F (X,X ;t), (16) 

where G F (X,X';t) = (2-rr)- 1 J dkexp[ik{X-X')-X(l- 
e- lka )t)] and the initial condition F(X', 0) = 5(X' - X ) 
has been noted. 

When t = 0, we have simply F(X, 0) =G F (X, X ; t = 
0) = 6(X — Xq). For large values of Xt, the function G F 
is dominated by small values of k; in this limit, using the 
saddle-point analysis, we find that Eq. (fl4l) is equivalent 
to the normal distribution: 



F(X,t\X ) 



1 



27T(T, 



exp 



(x - x - 

2a\ 



(17) 



where fi t = aXt, a t = ayXt, and F(X, t\Xo) stands for 
the distribution F(X, t) subject to the initial distribution 
F(X, 0) = S(X — Xq). The same result is obtained when 
the growth factor b and hence a are small; this allows us 
to expand the term exp(— ika) in the exponential inte- 
grand and to perform the gaussian integration. We may 
also expand G F directly as the series 

n (Y V'.-A f Ak(X-X')-Xt (MY* „-ikan 

G F (X ) X,t)=^j —e —e 



n>0 ' 



V- (A*) W -At 
2^ „! 6 



S(X - X' - an) 



= Y, F n(t)5(X-X' -an), 



(18) 



n>0 



which is nothing but the expansion in terms of the Pois- 
son distribution at given time. 

In consequence, the original size distribution assumes 
the log-normal distribution in the long-time limit: 



f(x,t\x ) 



1 



2ira t x 



exp 



(In a; — lnxo — (i t Y 



2a? 



(19) 

where f(x,t\xo) denotes the size distribution f(x,t) sub- 
ject to the initial distribution f(x, 0) = S(x — x$). Given 
an arbitrary initial distribution f(x, 0), the solution thus 
takes the form 



f{x,t)= dx'Gf(x,x';t)f(x',0) 



with 



Gf(x,x';t) = — G F (lnx, hix'; t) 



(20) 



(21) 
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Equation (|20|) . like Eq. (|4|, is invariant under the scale 
transformation x — » cx and x' — ?> cx'. 

We now consider the case r ^ 0. Expanding the second 
term on the right-hand side of Eq. (fTTj) : 



v 7 n— v 7 

(22) 

we write Eq. (fTTj) in the compact form 



Cf(x,t)=rg(x,t), (23) 
I 

f(x,t) = e~ rt [ dx' G f (x,x';t)f(x',0) + r 



where the linear operator L is given by 

d A °° (—1)™ / b \ n d n 

n— 

(24) 

Since the solution for r = 0, Eq. (fl"9j) . reduces to the 
Dirac delta function 8{x — xq) in the limit t — > 0, Green's 
function G(x,x';t,t') = G(x,x';t—t / ) for £ can be writ- 
ten in the form 

G (x, x'- t-t') = G f (x, x'; t-t') e-^-^e (t - t') , (25) 

with the Heaviside step function Q[t — t'). 

This leads to the general solution of Eq. (fTT|) in the 
form: 



* dt' e- T{t -^ [ dx' G f (x, x'; t - t')g(x', t'), (26) 



where Green's function reduces to the log-normal distri- 
bution in the long-time limit (t —¥ oo): 



creases exponentially in time; the second term thus be- 
comes dominant in the long time. 



Gf(x, x'; t) 



1 



2ira t x 



exp 



(lnx — lnx' — n t Y 



2a? 



The first term on the right-hand side of Eq. (121)1) , which 
describes the information on the initial distribution, de- 



A. Uniform-size production 

In case that new elements are produced in uniform size 
xq, i.e., g(x,t) — 5{x — xq), the asymptotic solution in 
the long-time limit reads 



f s (x) = f(x, *-*») = lim r / dt' e- r( *- t ' ) G/(a;, x ; t - t') = - / 



dk 



D ik ln(x/xo) 



x J 2n r + A(l - e- ika ) ' 



(27) 



I 

The integral over k can be performed with the help of be transformed into one over a unit circle in a complex 
the decomposition ka — 2im — where n is an integer plane: 



and < 9 < 2ir. Equation (|27l) then obtains the form of 
an infinite sum over integers n: 



2tt 



de 



-inO 



dz 



f — I_ ^ e 2i(TT/a)n \n(x/xp 



d6 e -( i / a ) 9ln ( x / x o) 
2tt r + A - Xe ie 



We notice that the integral over 9 is independent of n, 
and use the Dirac comb identity to write 



fsi x ) = — V<5(n-a l hx(x/x )) 
ax 

n 

= — 5 (n — a -1 ln(x/xo)) 



2tt 



l 2nr + X~ Xe %e ' J 2iirz r + A - Xz ' 

where we have set e %e = z. This integral vanishes when 
n is negative, i.e., for x < xq; otherwise Cauchy's residue 
theorem yields 

dz z- n _ A" 
2inz r + A - Xz ~ (r + A) n+1 
for n > 0. Setting finally x n = x e an = xq(1 + b) n , we 



d9 e"^/ a ) ln ( x / x o) 

2tt r + X — Xe 10 obtain the stationary distribution 



2n 



d9 



2nr + X- Xe 



if) ■ 



fs(x) 



r + X 



E - 

^ \ x Q 



7l>0 



5(x-x n ). (28) 



where it has been noted that the argument a 1 \a(x/x^) On the other hand, the stationary solution can be 

is constrained by each integer n. The integral over 9 can obtained simply by putting df/dt = and g(x,t) — 



5 



5(x — xq) in Eq. In this way, we still obtain the 

power-law distribution for x > x Q : 



the Dirac delta distribution, evolves to the power-law 
distribution in the uniform-size production process. 



fs(x) 



(29) 



B. Self-size production 



with the exponent a = 1 + ^rc^y > which should be 
exact in the stationary configuration. Equations (|28j) 
and (f2Uf give the same kind of power-law distribution 
in the long-time limit. It is indeed possible to recover 
from Eq. (|28[) the simple stationary power-law solution 
for small b or a. In such a continuous-growth limit, we 
also have r — > 0, with the ratio r/b and parameter A kept 
finite, since [ln(l+6)] _1 ln(l+r/A) and thus a would di- 
verge otherwise. In this case Eq. ([28)) may be approxi- 
mated by replacing the discrete sum over n by an integral 
over the variable u = an which becomes continuous in the 
limit a — > 0: 



ln(l + r/A) 

T 1 <r-^ / X \ ln ( 1 + b > 

/sW = — — r-/A — 5(an -\n(x/x )) 



r + \x ^ \xq 

n>0 v 



1-a 



T I f IX 

: / du \ — I S(u — ln(x/xn)) 

r + XaxJ a \x J 



XbxQ \xo 



(x - x ), 



When an element produces a new element of the same 
size as itself, the size distribution of newly produced el- 
ements is identical to that of the system, i.e., g(x.t) = 
f(x,t). In this case, the size distribution of newly pro- 
duced elements is not given externally, and Eq. (|26j) 
no longer provides the solution of Eq. (fTT|) but makes 
an integral equation for f(x,t). Remarkably, however, 
putting g(x, t) — f(x, i) in Eq. (fTTj) reduces the equation 
precisely to Eq. (j4]). Namely, the case of self-size produc- 
tion is effectively the same as the case of no production 
(r = 0), except for the increase of the total number of 
elements due to the production of new elements. In con- 
sequence, the size distribution of the system with self-size 
production is given by Eq. (|2"0|) . Starting from the initial 
configuration that all elements have the identical size xo, 
the size distribution develops to the log-normal distribu- 
tion: 



1 



2ir<j f x 



exp 



{\nx- n t y 

2a1 



(31) 



( 30 ) with fi t = ln(l- 



-b) Xt + Iiixq and a t = In (1+6) 
C. Binary fission 



Xt. 



where 9 is the Heaviside step function. This normalized 
power-law distribution has therefore the same exponent 
a = 1 + [ln(l+6)] _1 ln(l+r/A) w 1 + r/bX as the sta- 
tionary one given by Eq. (|29[) in the continuous-growth 
limit. Accordingly, our general solution given by Eq. 
(|2l))) demonstrates that the distribution, starting from 

I 



The binary fission by means of division in size may be 
interpreted as follows: The size of an element is reduced 
by the multiplicative factor /3 < 1 at rate r and a new 
element with the size reduced by the multiplicative factor 
1 — (3 is produced at the same rate r. As a result, the 
corresponding evolution equation obtains the form 



df(x,t) 
dt 



-(2r + X)f(x,t) 



—f(—t 
l + b J \l+b' 



r fx 

Ar 



i~p J Vi-/?' 



(32) 



r 



which describes both growth and reduction of the size 
due to binary fission. Although it can be easily ex- 
tended further to describe more complicated processes, 
we consider only one each kind of size-growth process 
and size-reduction process, which suffices for our study 

I 



here. Specifically, we suppose that the size of an ele- 
ment remains constant while the element gives birth to 
an offspring of the size proportional to its size with the 
multiplicative factor 1 — /?: 



df(x,t) 
dt 



-(r + X)f(x,t) 



1+b' 



,t + 



(33) 



with b > and < /3 < 1. Note that, in the half-size identical, 
division process (/3 = 1/2), Eqs. (1321) and (l33l) become 
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In contrast to the cases of uniform-size production and 
self-size production, the size distribution of newly pro- 
duced elements depends on that of the system and does 
not cancel out, respectively. Making use of the scaling 
symmetry of Eq. (|33l) . one is tempted to take a log- 
normal or a power-law distribution as an ansatz solu- 
tion. However, neither the log-normal distribution nor 
the power-law distribution can make a stable solution, 
as manifested easily by substitution. This may be in- 
ferred from the observation that the linear operator C 
for Eq. ([331) cannot be cast into the same form as that 
for Eq. (|24l) . Therefore we introduce another distribution 
observing the scaling symmetry, namely, the Weibull dis- 
tribution: 



with the shape parameter 7 and the scale parameter 77. 
There exist other distribution functions with the same 
symmetry, for example, the generalized exponential dis- 
tribution and the gamma distribution; they also fail to 
provide asymptotic solutions of Eq. (j33|) . 

Henceforth we denote y = \n(x/rj), u = 1 — (1 + 6)~ 7 , 
and v = 1— (1 — /3)~ 7 for convenience. Expanding (x/rj) 1 
at 7 = 0, 



'y 00 



(35) 



fe=0 



we write the time derivative of the Weibull distribution 
in the form 



15/ 7 / . 7 . 



OO 1 



(36) 







6) 7exp 


[-GDI 







while the right-hand side of Eq. 
(34) becomes 



33]) divided by f(x,t) 



-(rhs) = -(A + r) + A(l - «)e° + r(l - u)e u + A(l - «)e" ]T ^B k (v)( 7 y) k + r(l - u^c u £ -B fc ( u )( 7 y) fc (37) 
•' fe=i ' fe=i 



with the Bell polynomial given by (32j 

fl fc (z) = e-^^». (38) 

n=0 

Comparison of Eqs. (f5B|) and (|37p gives, for small 7, 

7 = {A [ln(l+6)] 2 + r [ln(l-/3)] 2 } + 0( 7 4 ) (39) 

to the zeroth order in y, and 

77 = 77 [A ln(l+6) + r ln(l-/3)] + 0( 7 ) , (40) 

to the nth order in y for all n > 1, where the expansion 
of v about 7 = 0, 

--1- (i + ft) =-Eir^« fc . ( 41 ) 
v 7 fe=i 

and similar expansion of u (with — f3 in the place of b) 
have been used. To the leading order, the solutions of 
Eqs. dSnj) an d P0)) under the initial conditions 7(0) = 7 
and 77(0) = ?7q read 

7 (<) = 7 ° (42a) 

"W=%e Dt , (42b) 

where C = A [ln(l+6)] 2 + r[ln(l-^)] 2 and D = 
Aln(l+6) + rln(l— j3). For this to be valid, one should 



I 

have sufficiently small 7 and y (or large 77) . In the asymp- 
totic limit the conditions of small 7 and large 77 require 
both C and D to be positive. Note that while C is always 
positive, D is positive only for r/X < — ln(l+6)/lri(l— 0). 
In conclusion, for r/X < — ln(l+6) / ln(l— 0), the Weibull 
distribution in Eq. ([34|) is asymptotically exact, where 
the shape parameter 7 decreases algebraically in time and 
the scale parameter 77 grows exponentially in time. 

D. Autocorrelations 

In the case of an asymmetric distribution of variables, 
there can exist power-law correlations in the variables for 
suitably adjusted parameters [33j. To probe the possibil- 
ity of long-time correlations in the evolution governed 
by Eq. (fTTj) . we consider the autocorrelation function at 
stationarity: 

_ ( X (T)x(T+t))-(x(T))(x(T+t)) 

vv(r)) - WT) 2 vV(r+t)> - W+W 2 

(43) 

where (x(t)) = Jdxf(x,t)x and (x(T)x(T+t)) = 
J dx dx' G(x,x' ;t)f(x' ,t)xx' with appropriate Green's 
function G(x, x'\ t) and r being large enough for the sys- 
tem to be in the stationary state. 

It is straightforward to compute Eq. (|43|) . which turns 
out to exhibit exponentially decaying behavior in the 
long-time limit. Specifically, in the case of the self-size 
production or no production (r = 0), Green's function 
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reduces simply to Eq. (|2Tj) : G(x,x';t) — Gf(x,x';t), and the resulting autocorrelation function reads 

C(t)~ex V (-^t\ (44) 
For the branching process, Green's function is given by 

I 

G(x,x';t) =e- rt G f (x,x';t)+r f dt'e'^-^Gfix^l-^x'^-t'), (45) 

Jo 



which leads to 



b^X+^r 
C{t) ~ exp 1 



(46) 



Finally, when there is a source producing new elements 
of uniform size, Green's function takes the form 



G{x,x';t) = e- rt G f (x,x';t) + r / dt'e'^'^Gfix, x ; t-t') 



(47) 



where xq is the size of newly produced elements. In this 
case the autocorrelation function depends on the produc- 
tion rate r: For r < 2bX + b 2 X, it reads 

C(t)~exp(-^*±t) (48) 

whereas for r > 2bX + b 2 X, we have 

C{t) ~exp[(-r + 6A)t]. (49) 

In all cases, the autocorrelations decay exponentially 
in time. We have also carried out numerical simulations 
to compute autocorrelation functions, to obtain results 
in full agreement with the analytical results listed above. 
It is thus concluded that long-time autocorrelations do 
not exist in the growth process considered here. This 
apparently implies that skew distributions do not neces- 
sarily result from long-time autocorrelations of the data 
variables. 



IV. NUMERICAL SIMULATIONS 

We carry out numerical simulations to check the an- 
alytical results in the preceding section. The master 
equation ([TJ is integrated directly to give the probabil- 
ity distribution for no-production, uniform-size produc- 
tion, self-size production, and half-size division processes. 
At each time step given by an integer multiple of At in 
simulations, a trial move is attempted according to the 
growth probability A = XAt and the production proba- 
bility R = rAt. To access the continuous-time case, At 
should be taken to be sufficiently small. Here we choose 




10 20 30 40 



x 

FIG. 1. Evolution of the size distribution in the case of 
no production, r = 0. Circles and squares plot simula- 
tion data for A = 0.2 and b = 0.025 at time t = 400 and 
500, respectively. Solid lines represent the corresponding 
analytic solutions, which are log-normal distributions with 
(/i,cr) = (1.975,0.2209) and (2.469,0.2469), respectively. The 
inset shows a (circles) and n (squares) versus time t, disclos- 
ing excellent agreement between simulation results (symbols) 
and analytical ones (solid lines). Error bars have been esti- 
mated by the standard deviations of the simulation data, and 
are smaller than the symbols. 



At = 0.001 and perform simulations of 10 2 different sam- 
ples, over which obtained results are averaged. Initially, 
all elements in each sample are taken to be of the unit 
size (xq — 1) and the number of elements is usually ad- 
justed so that each sample in the final data consists of 
N = 10 4 to 10 6 elements. These simulation conditions 
have been varied, to give no appreciable differences. 
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FIG. 2. Evolution of the size distribution in the case of 
uniform-size production. Circles and squares plot simula- 
tion data for (r, A, b) = (0.02,0.4,0.025) at time t = 500 
and (0.01,0.2,0.05) at t — 700, respectively, in the log-log 
scale. Solid lines represent stationary solutions, with the slope 
-a = -2.976 and -2, respectively, given by Eq. ([29]). The 
inset shows the exponent a of the distribution versus various 
parameters: Circles display a vs. the production rate r for 
A = 0.1 and b = 0.025; squares present a vs. the growth rate 
A for r = 0.1 and b = 0.025, and crosses a vs. the growth 
factor b for r = A = 0.1. Solid lines represent analytical solu- 
tions. Error bars estimated by the standard deviations of the 
simulation data are smaller than the symbols. 



When no new elements are produced, the simulation 
results are shown in Fig. [TJ and compared with the ana- 
lytic solutions. It is observed that the distributions ob- 
tained from simulations fit well to the log-normal dis- 
tributions described by the analytical expression in Eq. 
(|2U)) . The asymptotic solution also gives excellent es- 
timate for the time evolution of the mean and stan- 
dard deviation parameters (see the inset). This supports 
strongly the validity of the analysis in Sec. Ill, including 
the log-normal kernel in Eq. (I19|) and the asymptotic 
general solution in Eq. (|26l) . 

Figure [5] shows the simulation results for the case of 
uniform-size production, where the initial number of el- 
ements has been adjusted in such a way that the system 
ends with N = 10 6 elements for both cases, (r, A, b) = 
(0.02,0.4,0.025) and (0.01,0.2,0.05). For comparison, 
analytical results are shown together, manifesting power- 
law behavior. In particular, the inset exhibits the expo- 
nent a of the power-law distribution versus r, A, and b; 
there simulation results are compared with the analyti- 
cal solution given by Eq. (|28|) or Eq. (|29|) . As expected, 
excellent agreement is observed. 

For the self-size production process g(x,t) = f(x,t), it 
is expected that simulations give the same results as the 
no-production cases. To confirm this, we carry out sim- 
ulations of the system with the production rate r = 0.01 
and other parameters the same as those in Fig. Q] (i.e., 
A = 0.2 and b = 0.025). The resulting distributions at 
t = 400 and 500 (not shown) indeed turn out identical 
to those in Fig. [TJ This also implies that the contribu- 
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FIG. 3. Evolution of the size distribution in the half-size di- 
vision process with r — 0.002, A = 0.2, and 6 = 0.05. Circles 
and squares plot simulation data at time t = 200 and 300, 
respectively. Solid lines represent the fitted Weibull distribu- 
tion with (7, n) = (2.380, 6.834) for t = 200 and (1.842, 16.35) 
for t = 300, respectively. Error bars estimated by the stan- 
dard deviations are not larger than the symbols. The inset 
shows 7 (circles) and 77 (squares) versus time t, disclosing good 
agreement between simulation results (symbols) and analyti- 
cal ones (solid lines). 



tions of the size distribution of newly produced elements 
to the change in the size distribution of the system are 
described correctly by Eq. (fTT|) . 

Lastly, branching processes are simulated and the re- 
sults for the half-size division process (/3 = 1/2) are il- 
lustrated in Fig. [3] The resultant distributions fit well to 
the Weibull distribution. The shape parameter 7 and the 
scale parameter rj evolve in time, the agreement of which 
with the analytical results gets better as time goes on. 
This suggests that the branching process provides inher- 
ent mechanism of the Weibull distribution, the general 
origin of which has remained unresolved. 



V. DISCUSSION 

We have established a model for describing the size 
distribution in evolving systems and probed its solutions 
in several cases both analytically and numerically. Ob- 
tained are skew distributions including power-law, log- 
normal, and Weibull distributions, depending on the pro- 
duction rate, growth rate, and growth factor. Specifi- 
cally, power-law distributions are obtained when new el- 
ements of uniform size are produced; log-normal distribu- 
tions emerge in case that each element tends to produce 
either a new one of the same size as itself or no element 
at all. Further, branching processes give rise to Weibull 
distributions, suggestive of the origin of the Weibull dis- 
tribution. 

Emergence of the Weibull distribution in our model 
needs careful consideration: As the asymptotic expan- 
sion does not guarantee uniqueness of the solution, there 
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may exist other solutions which are distinguishable from 
our solution by the asymptotic expansion to the leading 
order. Further, there still lacks understanding of the time 
evolution of the Weibull distribution, e.g., ranges of the 
production and growth rates in which the Weibull dis- 
tribution provides quantitatively better fitting have not 
searched enough. The predicted time evolution of the 
shape and scale parameters is not exact, in contrast with 
the log-normal distribution for the self-size production 
process. Nevertheless, the physical meaning as to the 
emergence of the Weibull distribution is suggested con- 
vincingly and the accuracy of the asymptotic expansion is 
also supported to some degree by numerical simulations. 

With the insight gained from this study, we may un- 
derstand skew distributions observed in various systems. 
For instance, the power-law distribution of the popula- 
tions in the U.S. cities [34[ can be understood as an evolv- 
ing system, where small towns are created and then grow 
to larger cities. The number of passengers in an intra- 
urban subway system [35[ may also be treated by this 
model: The numbers of passengers making trips between 
two subway stations or using a subway station are indeed 
observed to follow log-normal or Weibull distributions. 
The size distribution of beta cells in pancreatic islets, 



which has been experimentally measured and theoreti- 
cally predicted [llj , is also described by a log- normal or 
Weibull distribution. It is expected that there are many 
additional evolving systems with skew distributions, ex- 
plained properly by this model. 

It is remarkable that ubiquitous emergence of such rep- 
resentative skew distributions as power-law, log-normal, 
and Weibull distributions is successfully described within 
a single framework. Applications to real systems exhibit- 
ing these skew distributions and more detailed under- 
standing of the emergence of the Weibull distribution are 
left for further study. 



ACKNOWLEDGMENTS 

One of us (MYC) thanks Institut Jean Lamour, 
Departement de Physique de la Matiere et des Materiaux 
at Universite Henri Poincare, where part of this work was 
carried out, for hospitality during his stay. This work was 
supported by the National Research Foundation of Korea 
through the Basic Science Research Program (Grant No. 
2009-0080791). 



[1] F. Lucchin and S. Matarrese, Phys. Rev. D, 32, 1316 
(1985). 

[2] M. Levy and S. Solomon, Physica A, 242, 90 (1997). 

[3] X. Gabaix, P. Gopikrishnan, V. Plerou, and H. E. Stan- 
ley, Nature, 423, 267 (2003). 

[4] G. K. Zipf, Human Behavior and the Principle of Least 
Effort (Addison- Wesley, Cambridge, 1949). 

[5] V. Pareto, Cours d'Economie Plitique (Droz, Geneva, 
1896). 

[6] G. U. Yule, Phil. Trans. Roy. Soc. London. Ser. B, 213, 
21 (1925). 

[7] R. Gibrat, Lea Inegalites Economiques (Sirey, Paris, 
1931). 

[8] A.-L. Barabasi and R. Albert, Science, 286, 509 (1999). 
[9] W. W. Fleming, D. P. Westfall, I. S. de la Lande, and 
L. B. Jellett, J. Pharmacol. Exp. Ther., 181, 339 (1972). 
[10] K. Oexle, J. Theor. Biol., 190, 369 (1998). 
[11] J. Jo, M. Y. Choi, and D.-S. Koh, Biophys. J., 93, 2655 
(2007). 

[12] E. Limpert, W. A. Stahel, and M. Abbt, BioScience, 51, 
341 (2001). 

[13] F. Black and M. Scholes, J. Political Econ., 81, 637 
(1973). 

[14] W. J. Reed, Physica A, 319, 469 (2003). 
[15] M. Mitzenmacher, Internet Math., 1, 226 (2004). 
[16] W. Weibull, ASME J. Appl. Mech., 18, 293 (1951). 
[17] T. W. Peterson and M. V. Scotto, Powder Tech., 45, 87 
(1985). 

[18] G. R. McDowell, M. D. Bolton, and D. Robertson, J. 

Mech. Phys. Solids, 44, 2079 (1996). 
[19] V. Papadopoulou, K. Kosmidis, M. Vlachou, and 

P. Macheras, Int. J. Pharmaceutics, 309, 44 (2006). 



[20] J. Beirlant, Y. Goegebeur, J. Segers, and J. Teugels, 
Statistics of Extremes: Theory and Applications (Wiley, 
Chichester, 2004). 

[21] E. Bertin and M. Clusel, J. Phys. A, 39, 7607 (2006). 

[22] K. Kosmidis, P. Argyrakis, and P. Macheras, Pharma- 
ceutical Res., 20, 988 (2003). 

[23] K. Kosmidis, P. Argyrakis, and P. Macheras, J. Chem. 
Phys., 119, 6373 (2003). 

[24] J. Jo, J.-Y. Fortin, and M. Y. Choi, Unpublished. 

[25] W. K. Brown and K. H. Wohletz, J. Appl. Phys., 78, 
2758 (1995). 

[26] K. Ohashi, L. A. N. Amaral, B. H. Natelson, and Y. Ya- 

mamoto, Phys. Rev. E, 68, 065204 (2003). 
[27] C. Braun, P. Kowallik, A. Freking, D. Hadeler, K.-D. 

Knifiki, and M. Meesmann, Am. J. Physiol. Heart Circ. 

Physiol., 275, H1577 (1998). 
[28] C. L. Ehlers, J. Havstad, D. Prichard, and J. Theiler, J. 

Neurosci., 18, 7474 (1998). 
[29] M. Sammon and F. Curley, J. Appl. Physiol., 83, 975 

(1997). 

[30] M. Y. Choi, H. Choi, J.-Y. Fortin, and J. Choi, Euro- 

phys. Lett., 85, 30006 (2009). 
[31] H. Risken, The Fokker-Planck Equation: Methods of 

Solution and Application (Springer- Verlag, New York, 

1996). 

[32] S. Roman, The Umbral Calculus (Academic Press, 
Florida, 1984). 

[33] B. Podobnik, P. C. Ivanov, V. Jazbinsek, Z. Trontelj, 
H. E. Stanley, and I. Crosse, Phys. Rev. E, 71, 025104 
(2005). 

[34] M. E. J. Newman, Contemp. Phys., 46, 323 (2005). 
[35] K. Lee, W.-S. Jung, J. S. Park, and M. Y. Choi, Physica 
A, 387, 6231 (2008). 



